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ABSTRACT . 


This work examines the three-dimensional dynamic response 
of a nonlinear fast reactor with temperature-dependent feed- 
back and delayed neutrons when subjected to uniform and local 
disturbances. The finite element method was employed to re- 
duce the partial differential reactor equation to a system of 
ordinary differential equations which can be numerically in- 
tegrated. A program for the numerical solution of large 
Sparse systems of stiff differential equations developed by 
Franke and based on Gear's method solved the reduced neutron 
dynamics equation. Although a study of convergence by refin- 
ing element mesh sizes was not carried out, the crude fi- 
nite element mesh utilized yielded the correct trend of neu- 


tron dynamic behavior. 
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I. INTRODUCTION 


The nuclear reactor with delayed neutrons and tempera- 
ture-dependent ae at is a nonlinear system, whose response 
to uniform and local disturbances differs greatly from that 
Of a linear reactor. The prompt temperature feedback model 
and one average group of delayed neutrons were incorporated 
in this work. Practically all neutron dynamics analysis to- 
day deals with two-dimensional geometry, implying a symmetry 
in neutron dynamics behavior during transient [1,2]. This 
Symmetry assumption, however, could be unrealistic in the 
Safety analysis of nuclear reactors. A unique feature of 
this work is the consideration of the three-dimensional depen- 
dence of the neutron flux. No symmetry assumptions were im- 
posed on the problem. This work investigated neutron dynamics 
under uniform and local disturbances in the core. A uniform 
initial flux throughout the interior of the reactor was im- 
posed. The finite element method (FEM) was employed to solve 
this nonlinear initial-boundary value problem. The FEM is 
quite effective in handling discontinuous forcing functions 
thereby making it particularly suited for examining the ef- 


fects of localized perturbations and space-dependent feedbacks. 
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II. THE NUCLEAR REACTOR WITH TEMPERATURE 
DEPENDENT FEEDBACK AND DELAYED NEUTRONS 


Boy THE PHYSICAL. SYSTEM 

The system under consideration is a fast reactor of cy- 
lindrical geometry that is composed of two different regions. 
The reactor core or region I is cylinderical in shape and is 
fueled by U-235. Region II or the reflector region complete- 
ly surrounds the core and is composed of U-238. Both regions 
were assumed to be homogeneous. Table I lists the physical 
properties and geometry for each region. A schematic of the 
fast reactor geomenry is shown in figure 1. 

Temperature and delayed neutron effects were taken into 
account. Also, a one-velocity Or one-group model was assumed, 
thereby making the velocity independent of spatial or tem- 
poral effects. The delayed neutrons were considered only in 
region I since region II was assumed to be a non-multiplying 
medium. The temperature effects were also assumed to be only 
in the core region. 


In general form, the one-velocity neutron diffusion equa- 


tion is [3] 

aN( r,t) _ Si ie - 

= vDV?N - Z.vN + S(r,t) (1) 
where ae Ce | 


D = neutron diffusion coefficient 
=Z. = macroscopic absorption cross-section 


N(r,t)dV = number of neutrons in a volume element dV at 
a point rat time < 
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SYMBOL 


TABLE I 


Physical Constants 


DEFINITION 
radius of Region I 
total reactor radius 


height of Region I 
total reactor height 


neutron speed 


core neutron diffusion coefficient 


reflector neutron diffusion 
coefficient 


core neutron absorption cross 
section 


reflector neutron absorption 
cross section 


number of neutrons per fission 
critical fission cross section 


delayed neutron fraction; 
dollar reactivity 


fission energy 

modified convection heat 
transfer coefficient 
temperature coefficient 


abundance-weighted mean decay 
constant 


12 


VALUE 


60 cm 


90 cm 


160 cm 
220 cm 


u.sx10/ cm/sec 
0.913 cm 
1.200: cm = 


0.01401 om 


0.0040 om 


2.904 


0.005736 an + 


0.00642 


7.652x10 
eal/fission 


0.0632 


cal/(Cem- 


sec °c) 


-~0.004/°C 


0.4350 sec + 


4 to \) 
# =P he : ae a ' 
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ual 7 a #7, 
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a. top view 


fan ~ Core b. front view 
II - reflector 


Figure 1. Schematic of cylindrical resgeror. 
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vDV2NdV_ 


number of neutrons diffusing into dV per 
unit time at time t 


rVN dV = number of neutrons absorbed in dV per 
unit time at time t 
S(r,t)dV = number of neutrons produced in dV per 


unit time at time t 


The neutron number density is related to the neutron flux by 


the expression [4] 

Sie,t) = eh r.t) (2) 
where (r,t) is the flux at time t . The neutron diffu- 
sion equation in terms of the flux is depicted by 


. oett st) = DV*o - £.9 + S(¥,t) (3) 


B. PROMPT AND DELAYED NEUTRONS 

The source or production term in equation (3) is composed 
of the contributions of the prompt and delayed neutrons. The 
majority of the fission neutrons are prompt neutrons that 
appear almost instantaneously (within ae second) on fission. 
Assuming a fast neutron non-leakage probability of unity, the 


prompt neutron source is described by 


S(t) = (TEST Kites t)eOlreer (4) 
where 
RoLh,t) = infinite medium multiplication factor 
8 = total fraction of delayed neutrons 
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The macteus of delayed neutrons is very small (note 
that 8 = 0.00642 from Table I). However, they have a very 
Significant effect on the reactivity because their mean life- 
times are long. Without delayed neutrons, reactor control 
would not be setad DUES. These delayed neutrons are born in 
the decay by neutron emission of nuclei produced following 
the B-decay of certain fission fragments. For example, the 


87 leads to Kreo 


B-decay of the fission fragment Br plus a 
neutron. Nuclei such as pre? whose production in fission 
eventually leads to the emission of a delayed neutron are 
called delayed neutron precursors [4]. 

There are six main groups of delayed neutrons. Each 
group is classified according to its decay constant. The de- 


layed neutron source term is portrayed by 


Sy(r.t) = Ci(r,t) a, (5) 
i=] 
where 
th 
d; = decay constant of the i group 
C.(r,t) = density of the jth precursor 


Assuming that the fission fragments do not migrate appre- 
ciable distances and assuming a non-circulating fuel reactor 
[3], the precursor density is delineated by 
rR = — 

= B,K,2,¢ - AC, (6) 


where B is the fraction of delayed neutrons of the yen 
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group. The solution to the precursor equation is in terms 


of a time integral expressed by 


Gir.t) = et, f a hgh ine Tie ce eherece iat! (7) 


0 
Inserting equation (7) in equation (5) yields the delayed 


neutron production term as 


° 3 (8) 
For convenience the six main groups of delayed neutrons 
were considered as one group. This was accomplished by us- 


ing the abundance-weighted mean decay constant defined by 


= 1 6 
Se (9) 
i=1 


This is a reasonable approximation since many of the impor- 
tant phenomena of nuclear reactor dynamics can be character- 
ized satisfactorily by combining all the emitters or precur- 
sors into one, two, or, at most, three effective groups [3]. 
Replacing dj with the abundance-weighted mean decay con- 
stant showed the delayed neutron source term to be 


ae 
pre) ee, fet kK (Ft) o(F,t) dt! (10) 
O 


C. DOPPLER EFFECT AND TEMPERATURE FEEDBACK MODEL 
The Doppler effect in fast reactors is due to the tempera- 
ture broadening of many closely spaced high-energy resonances 


in both the fission and parasitic-absorption cross-sections. 
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These resonances basically mean that as the temperature in- 
creases, the number of neutrons that are absorbed increases 
and, similarly, the number of fission events increases. 
These nonproductive and productive processes compete in a 
complicated manner, and the net effect may be either an in- 
crease or decrease in reactivity [3]. For most fast reac- 
tors, the effect is negative, and it will be so assumed in 
this work. 

The reactivity change with respect to the fuel tempera- 


ture change is modeled by [2] 


where a,b, and care parameters determined from experimen- 
tal or neutronic calculations, m is an integer, and 7 15s 
the current fuel temperature. For most fast reactor cores 
(ceramic-fueled cores), TdK /dT is very close to being con- 
stant over a wide range of temperatures. Therefore, it is 
assumed that a and care zero and the Doppler coefficient 


is expressed by 


The following initial conditions were used: 


we Beas (12a) 
T (gee Pee (12b) 
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The reactivity model is then expressed by 


K(P,t) = Ko + b gn (=) ae i: 
fe) 
where 
K2 = multiplication factor at steady-state 
v = number of neutrons produced per fission 
Fe = original fuel temperatures 


= Doppler constant 


The definition of Ke 1's 


5.” 
Vv 
adbeast (14) 
a 
where 
* 
Le = critical fission cross-section 
= macroscopic absorption cross-section of the 


core 


The rise in fuel temperature is depicted by 


This is also described by the integral [5] 


3 
(r,t) = f(t-t') w(r,t)dt! (16) 
0 


where f(t-t') is the feedback kernel and is dependent upon 


the type of temperature feedback model used. w(r,t) is the 
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rise in neutron flux above the steady state and is expressed 


by 


where o'r) is the neutron flux at steady state. 


There are three types of temperature feedback models that 
can be considered here. The first is known as "Newton's 
feedback" which determines the reactor temperature by Newton's 
law of cooling. The second temperature feedback model is 
called the adiabatic feedback model. This represents the 
temperature for the loss of coolant case. The third model is 
the prompt temperature feedback model in which the fuel tem- 
perature follows the behavior of the neutron flux without de- 
lay [5,6]. The prompt feedback model was employed in this 
analysis. The feedback kernel for this temperature feedback 


model is 
i ea ' 
f{t=t! j= a 6lt=£ (18) 


where K is an energy production operator with units of °C 
per unit flux and y , a dimensionless quantity, is related 
to the mean time for heat transfer to coolant. 

Inserting equation (18) in equation (16) and performing 


the integration produced a rise in temperature of 
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From equation (19) the temperature of the fuel is 
ee a - 
T(#,t) = Sy(tt) + 1,(F) (20) 


The ratio of the current temperature to the steady-state tem- 


perature is therefore 


| Mee; Pees 1 (ona 


In this work te was considered to be constant throughout 
the core. Incorporating equation (21) into equation (13) 


gave the reactivity model of 


K(7,t) = Ko +b en [— y(F,t) + 1] (22) 
Y'o 


D. FIELD EQUATIONS 

Before establishing the field equations it is desirable 
to express equation (3) in terms of the rise in flux. Using 
equation (17) in equation (3) and grouping terms yielded the 
following diffusion equation: 


t 
ow _ ha be t Ha x 
ae 7 COV y-t,y +(1-8)K 2b + BAZ, e 
O 


Mt-t dy vate] 


<|- 


t 


+ LDV O57 Toot (1-8) KF 64,7 an, f e 


-A(t-t' kK odt'] 


(23) 
The second bracketed term in equation (23) is identically 


equal to zero since it is the steady state portion of the 
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diffusion equation. The rise in neutron flux above its 


steady state value is therefore expressed, for the core, by 
yap = OVy - Ely + (1-8) KE.y 
+ BAL, fp galt) Kv dt' (24) 
and, for the reflector region, by 
1 3Y = DV7y - Ey (25) 


Inserting the reactivity model into equation (24) yielded 


ae 2 = = ° 
oa DV-p zo + (1-8) Za Key 


t 
+(1-8)bz, cnt a et goatee fg Ee nae 
ie) 
= E -k(t-t') Ky 
+ Biz.b ie e [anor + 1)] vdt'} (26) 


6) 


For the reflector, the last four terms of equation (26) are 
zero. The effects of the temperature on the delayed neutrons 
were neglected in this work. The field equations can now be 
expressed, for the core, as 

oy _ 2 z ; 

oe 7 ON + [vz. v(1-8) vig]y 
+ [-(1-8)ve_b] Len(S@ + 1)]y 

a rae 


+ [-pivo_v] hie (teEC dy eee (27) 
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and, for the reflector, as 
P) 
+ - vDV7y + vE.y = 0 (28) 


The non-linear terms of the core field equation will be 
linearized accordingly. In more compact form, equations (27) 


and (28) became 


t y 1 
ae = clV*» + c2y + c4fun (St +1) wrest s geen dydt' J=0 
(@) 0 
(29) 
and 
3 
= - clv2w + c2y = 0 (30) 


where the Meanings of the coefficients cl, €2, etc... are 
obvious for the core and reflector. 

Equations (29) and (30) were subjected to the following 
conditions: 

boundary condition: 


Wa(F-t) = 0 (31) 


where rp are coordinates of points on the outer surface of 
the reflector and the subscript R_ refers to the reflector 
eon. Loy. of .f lux. 


vp(r.t) = vi (ry.t) (32) 


where ry are coordinates of points on the core-reflector 


interface and the subscript c_ refers to the core. 
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Tit. FINITE ELEMENT 


A. INTRODUCTION 

The application of the finite elements of nonlinear con- 
tinua has been mostly in the field of solid mechanics. 

Prior work [1] using the finite elements of nonlinear con- 
tinua on a nonlinear reactor dynamics problem has been suc- 
cessful. The FEM in this work was utilized to reduce a non- 
linear partial differential equation of the nuclear reactor 
to a system of nonlinear ordinary differential equations in 
time. The time integration was accomplished by using a com- 
puter program for the numerical solution of stiff differen- 
tial equations developed by Franke [7]. In order to minimize 
computer storage requirements, an optimum compacting scheme 
(OCS), described by Ref. 8, was adopted. 

The finite element models of operator equations are gen- 
erally classified into three categories: a) the variational 
finite element models such as the Ritz method, b) the weigh- 
ted residuals method such as the method of Galerkin, and c) 
the direct finite element models which are not based on 
functional minimization. From experience in structural 
mechanics, the most effective method for generating accept- 
able finite element models of nonlinear equations is the 
Galerkin method [1]. This work adopted the method of Galer- 
kin in a finite element approximation over the spatial 


domain of the field equations. 
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B. THE METHOD OF GALERKIN 

The Galerkin method is a special case of the method of 
weighted residuals. It involves a rational choice of 
weighting function that is consistent with the type of fi- 
nite element approximation considered. Indeed, the weight- 
ing functions chosen are the basis or shape functions em- 
ployed in the finite element approximation. There are two 
favorable characteristics of the Galerkin method which 
makes it attractive. The first attribute is its amenability 
to integration by parts. This supplied the freedom of using 
a lower order finite element than might be otherwise possible. 
The second favorable characteristic of the Galerkin method 
is that the symmetric operators in the field equations trans- 
form into symmetric matrix operators. Both these attributes 
are attractive for computational purposes. 


Consider the initial-boundary-value problem 


ae = fy - f(F,t) (33) 


where ET contains the nonlinear operators. According to 
the spatial finite element discretization, the solution of 
equation (33) is in the form of an union of an N-term approx- 
imation given by 


A 
v(#,t) = v(F,t) = U Yt) BIL) y J=1,2,...,8 aD 


where N is the number of system degrees of freedom (i.e., 
number of coordinates), and Gj(r) are the system or global 


basis functions which span the space of the approximate 
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solution w(r,t) [1]. The Einstein summation is used. The 
global basic functions tte “pyramid functions", each of 
which has a eee er functional description over a sub- 
domain of the system and is zero elsewhere [9]. The unknown 
coordinate functions v(t) are the time-dependent magni- 
tudes of the approximated flux w(r,t) and/or its deriva- 
tives at discrete nodal points [10]. 


The residual function, R(r,t) , is defined such that it 


is identical to zero when w(r,t) is equal to the exact solu- 


tion. The residual function is expressed by 


Re > eee (35) 


The Galerkin orthogonality condition (using the basis func- 
tions as weight functions), when applied to the residual func- 


tion, requires that 


GR(r,t)dVol = Os T2296 oh on (36) 


Yo] 


From the field equations, the residual function for the core 


R(r,t)= Ae clv2y + c2u + cALen (st + 1)]v 


3 0 
t - 
+ est [ en (t-t') faery (37) 
ie) 
and for the reflector 
R(r.t) = 2 - clv2b + c2y (38) 
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Using equation (34) and applying Galerkin's orthogonality 


condition produced a core equation of 


| EXCNOREAGESE + Gy(c2-y), 


Vol 
K 
*Len( or Grby + bh) Gy(c4-y) | 
t... 
ef eA(t-t (5+), dt']G 5} = 0 (39) 
) 
where kaa eee | 


The reflector has a similar equation without the nonlinearity 


and is expressed by 


rei ce) - VG (cl+v) | ‘ Gy(c2-) 3 dVol = 0 (40) 
Vol 


C. THE ELEMENT 

A three-dimensional quadratic isoparametric element was 
employed in this work. The parent element is a triangular 
prism or solid wedge with straight sides. The element shape 
functions are expressed in terms of area coordinates in the 


plane of the triangle and by an isoparametric coordinate 
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along the prism axis. Figure 2 shows the parent element. 
This element was chosen because of the ease with which it 
fits the cylindrical structure when it is transformed into 

a curved element. This type of element has been used before 
as filler elements [11]. 

The area coordinates are defined by area ratios. Con- 
Sider the triangle shown in figure 3. An arbitrary point P 
within the triangle defines three subareas designated by A,» 
Ay s and A. The ratio of each of the subareas to the total 
area is known as an area coordinate. In equation form, the 


area coordinates are 


Ly = A,/A (41a) 

L, = An/A (41b) 

L3 = A/A (416) 
where A is total area of the triangle. Ly> Los and L 3 


are the natural coordinates for a triangle. The requirement 
that the sum of the subareas be equal to the total area is 


Obviously satisfied by the identity 


te eg eS ee (42) 


In the plane of the triangle, only two of the area coordi- 
nates are independent. 

The isoparametric coordinates are best visualized by con- 
sidering the rectangular prism shown in figure 4. Isopara- 
metric coordinates are normalized coordinates such that 
their values on the faces of the rectangle are + 1. The &nz 


axes are in general not orthogonal. They are orthogonal 
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Figure 2. Quadratic triangular prism parent element 
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Figure 3. Definition of area coordinates 
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Figure 4. Isoparametric coordinates 
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only in the special case of a rectangular prism element 
[12]. The element basis functions use the cz coordinate 
in the prism axis and the Lys Lo» and L3 coordinates in 
the plane of the triangle. 

The parent element, as shown in figure 2, has 15 local 
nodal points around its periphery. As such, there are 15 


basis functions which are given below [11]: 


forner nodes: (nodes 1, 3, 5, 10, 12, 14) 
Ny = 5 by (2Ly-1) (142) - F Ly (1-8?) 
Ng = x bp (2by-1) (142) = F Ly (1-2?) 
Ne = by(2l4-1)(1+c) - Lg (1-2?) 
No gt y by (2by-TW) (1-8) = F by (1-2?) 
Not gy bp (2by-1) (1-5) - 5 by (1-0?) 
Nigt g by(2bg-1) (1-5) - 5 bg (1-2?) 


Midside nodes of rectangles: (nodes 7, 8, 9) 


Midside nodes of triangles: (nodes 2, 4, 6, 11, 13, 15) 


No = 2L,L,(1+s) 
m= 2Lyl,(1+2) 
Ne = 2L,L,(1+c) 
Nj45 2L,L,(1-<) 
N1 37 2L5l4(1-) 
Nis 2L4L,(1-<) 
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The coordinates of each local: node, in terms of L a L 


Ro ae ee 
and ct are listed in table II. These element basis functions 
< N > define the geometry of the element. Note that they 


satisfy the relationship 


1, @t node 7 


0 , at other nodes (43) 


On the element level, the variation of the unknown func- 


tion w is approximated by 


ME fl te os 2 (44) 
where © < N'> = row vector of element shape functions 
{y}© = column vector of time-dependent nodal 


values of wv 


To satisfy continuity requirements, the shape functions 
< N'> have to be such that the continuity of the unknown 
function wy is preserved in the parent coordinates [11]. 
The shape functions < N > which characterize the element 
geometry and the shape functions < N'> which describe the 
unknown function do not necessarily have to be the same. 
There is no requirement that the nodal values be associated 
with the same nodes which were used to define the element 
geometry, though in practice it is often the case. Consider 
for example, the illustrations in figure 5, If the nodes 


defining the element geometry and the nodes defining the 
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TABLE II 


Coordinates of Local Nodal Points 


Local Node 
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ay 
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aw 
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(a) Isoparametric 


(b) Super-parametric 


(c) Sub-parametric 
@ - geometry nodes 


O - variable function nodes 


Figure 5. Element classification 
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unknown function are identical, the element is known as an 
isoparametric element. This means that the shape functions 
describing the geometry and the shape functions describing 


the unknown function w are equal, or 
< W's -= <> fi > (45) 


If there are more nodes defining the geometry than nodes de- 
fining the variable function, the element is called a super- 
parametric element. Using more nodal points to define the 
unknown function than to describe the geometry of the element 
results in a subparametric element [11]. This work utilized 


the isoparametric element classification. 


ie otVisSion OF THE SYSTEM INTO ELEMENTS 

In three-dimensional space, the division of the sys- 
tem into discrete elements is difficult to visualize. It is 
virtually impossible to show every nodal point of the system 
in one schematic. To present a clear view of the discretized 
domain, a "layer" approach was adopted. The first finite 
element grid or mesh employed here consisted of 128 elements. 
Under this grid (mesh I), the reactor was divided into four 
layers as shown in figure 6. Each layer was composed of 32 
elements. The first and fourth layers were each 30 cm in 
height and each contained entirely reflector elements. The 
second and third layers were each 80 cm in height and together 
they encompassed the entire core plus the remaining reflector 
elements. Each layer, in turn, was partitioned into three 


horizontal (xy) planes. The top plane included all the global 
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Figure 6. Layers of mesh I 
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or system nodes corresponding to local or element nodes 1 
through 6. The middle plane contained all the system nodes 
corresponding to the local nodes 7, 8, and 9. System nodes 
corresponding to element nodal points 10 through 15 comprised 
the bottom plane. To fix ideas, the three nodal planes of 
the first layer of mesh I are shown in figures 7, 8, and 9. 

In this work there was only one curved side which was in 
the plane of the triangle, as shown in figure 10. The first, 
seventh, and tenth element nodes were each arbitrarily assigned 
to have as an opposite side the curved side of the triangle. 
The remaining local nodes in each of the respective planes 
of the element were then numbered consecutively in the counter- 
clockwise direction. 

At this point it is appropriate to define connectivity. 

The connectivity of an element is a row vector array that 
relates the local nodes to system nodal points. The connec- 
tivity lists the system nodal points that are "connected" to 
form the element domain in the sequence of local nodal number- 
mag. thus, the connectivity matrix: for mesn [ is a matrix 

of size 128 x 15. To illustrate, the connectivity of element 


number 106 is 
mes5, 210, 193, 194, 1195,c2]1, 266, 176, 177, 263, Pe, 10, 11, tee 


A second finite element mesh (mesh II) consisting of 192 
elements was developed. It contained the same number of ele- 
ments per layer as mesh I. However, mesh II has six layers. 


The second and third layers of mesh I were each divided in 
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Figure 8. Middle nodal plane of the first layer of 
mesh I, z = 95 cm 


a3 


yeGryun tneme's - sedmult 


* 3 ,? @taa 


e weyet teh? ond mprerets taken Steprm .@ er vwoeta 
mo 


d 


Bu2 2B 282 


307 


306 


305 


304 


303 


298: 997 296 


Number - element number 


Figure 9. Bottom nodal plane of the first layer of 
mesh I, z = 80 cm 
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half, thus forming the two additional layers of mesh II. 
The connectivity matrix and the coordinates of each global 
node for mesh I are given in Appendix A. Appendix B lists 
the connectivity matrix and nodal coordinates of mesh II. 
The reader can construct the different layers and elements 
without much difficulty by using the nodal coordinates and 
the connectivity matrix. A mesh generator was not utilized 


in this work. 


E. COORDINATE TRANSFORMATION 

The use of a quadratic or higher order element permits 
the transformation or mapping of the straight-sided parent 
element into an element with curved sides. Distorted or 
curved elements provide a better fit to curved domains than 
linear elements, and thus a smaller number of elements is 
required to represent the structure adequately. 

The transformation from cartesian coordinates to curvi- 
linear coordinates can be accomplished by employing a one-to- 


one correspondence defined by [11] 


eB 


The element shape functions are utilized to achieve this 


(46) 
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transformation via the relation 


4 2 Nix, 
a 
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where n° is the number of element nodes, and the element 
shape functions N; are in terms of local coordinates. Each 
set of local coordinates corresponds to only one set of car- 
tesian coordinates. 

In performing the transformation, the vandal eit re- 
quirement must be met. The transformation into the new, 
curved elements should leave no gaps between adjacent elements. 
If two adjacent elements are generated from parents in which 
the element shape functions satisfy continuity requirements, 
then the curved elements will be contiguous [11]. For the 
isoparametric element, uniqueness of coordinates ensures com- 
patibility. Continuity is assured when adjacent elements 
are given the same sets of coordinates at common nodes. 

Since the element shape functions are in terms of local 
coordinates, an element of volume, dxdydz, must be transformed 
into an element of volume expressed in local coordinates. 
This-is achieved through the use of the Jacobian matrix de- 
fined below. Using the chain rule, the relationship between 


E,n,t and a corresponding set of cartesian coordinates 
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The Jacobian matrix [J] is defined as 


ne) 5yAs age 
de * 96 7 ge 
ff) ae. See a 
ax ay) * dy 
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From equation (47), the Jacobian matrix becomes 


JN. ON. ON. 
TM aera na eae 
ON. oN. ON. 
ss nies | ee call 
[J] : an Ay ; an Yi? : an ai (50) 
oN. oN. ON. 
eee tt Se goes Brees 


The determinant of the Jacobian is used to transform the vol- 
ume of element in cartesian coordinates to local coordinates. 


For a volume of element [11], 
dxdydz = det [J] d&dndz (51) 


Note that the determinant of the Jacobian is a variable for 
elements of curved geometry. Only in the case of straight- 
sided elements is the determinant of the Jacobian a constant. 
In the plane of the triangle, the area coordinates (Ly, 
L5>L3) number one more than the cartesian coordinates (x,y). 


Thus, L, is defined as a dependent variable. This establishes 
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the origin of the &n coordinate system at corner point 3, as 
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illustrated in figure 11. Recall that the En axes need 


not be orthogonal. As such 


— = Ly (52a) 
n= Lo (52b) 


Using equation (42), 


L4 =1-E&-7n (52c) 


Applying the chain rule yields 


oN. : ON. dL, ON. dL. : ON. dl, (53a) 
dg 3 1 3 P) 9 ag dL, era 
ON. : ON. dL, n ON. dL, ‘ ON. dL. (53b) 
an R) 1 an Ll, n dL, on 


Using equations (52a), (52b), and (52c) in equations (53a) 
and (53b) gives 


oe 
mL, aL, (54a) 
aN, aN. aN, 

mae st, 3, (54b) 
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Now the Jacobian matrix can be evaluated using the shape func- 
tions N. = Ni (L),l,5L,,¢) . Inserting equations (54a) and 
(54b) in equation (50) produces 


eb oe peleeels pads ae 
Ula = a) x 9 es ray Ear 
dL, = aL 4° aL; aL’ i aL, L4G 
i #52 3 jaar 4 9n, te 3 
ON oN. oN. 
aati pelt pie! 
(Fe )x. > (ss )y; > ot )z. 
(55) 


To summarize, suppose it is required to transform the 


integral 


oe ne Fi (bs host Sse), dydydz (56) 
Vol 


to an integral entirely in terms of local coordinates. The 


determinant of equation (55) is then utilized to give 


+1 1-t 
Top 
ey ff F'(L, Lo sb9.¢)det[J(L, LL, 2) JdL, dL de 
Fifa /o a} 


Equation (56) is now in a form suitable for numerical inte- 
gration. 

The development of coordinate transformation to this 
point has been under the general assumption that all the 


sides of the element are curved. In this work, the only 
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curved side is in the plane of the triangle. The sides of 
the element along the prism axis are straight. As such, a 
modified Jacobian matrix can be employed, thereby reducing 
the number of calculations to be performed. This modified 


Jacobian is the 2x2 matrix defined by 


oN ON. oN. ON. 
ae 1) x Fiat 3 lyy 
qeetqge F371” 8k, ett 
os ON. oN. ON. ON. 
[J ] = EGE, = ae ’ aes = aL, Yi (58) 


Along the prism axis or ¢ direction, it can be shown that 


rm|> 


dz = = dt (59) 


where h = height of the element. The volume relationship 


is then given by 


dxdydz = 3 det[J ] dL, dlc (60) 


The integral of equation (56) then assumes the form of 


ett 
1 
2 0 ; “ 
r-Sf ff F(L,sb5sb3.¢)det[J (L,.Ly.L3.c)]dL.dl,,de (61) 
a oO 'G 


Equation (61) is the basis of numerical integration applied 


in this work. 


F. CONSTRUCTION OF ELEMENT MATRICES 
The system matrix operators can be constructed through 


the use of the global basis functions G, or through the 
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application of the element shape functions. Although the 
global basis functions were used to demonstrate the method 
of Galerkin, the construction of the system matrices in this 
work was achieved through element considerations. The solu- 
tion of the unknown variable w within the element domain 


was approximated by 
~e 
igi Pal Pee tas a ce ee (62) 
i=] 


where n° is the number of element nodal points, N, are 

the element shape functions, and v5 (t) are the time-depen- 
dent nodal magnitudes of ve . The element contribution to 
the system matrix operators is defined by Galerkin's ortho- 


gonality condition expressed by 


és N, R° dVol = 0 , j=1,2,...n° (63) 


where R° is defined by replacing v with ye in equations 
(37) and (38), and the integration is over the element volume. 
Using equation (62), the element contribution is portrayed 
for the core by 
oe 2 ene owe 
ts NN, dVol]{b, ce)it J N.VN,dVo1]{ (cl +p ae NN;dVoT]{(c2-u") ;} 


K = e 
+ [ J N.N.2gn(1 + ——2N, wv, -(t))dVol] {(c4-v~).} 
J TAl ¥hy k k <k 1 


t x ' 
+L SN Ndvol] tf en A(t-t') (es .y(t") ).dt'} = 0 (64) 
Vol 0 
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mneres isjske= 1525.25 35955,° and cly e2gretey staré constants 
at node i. The bracketed expressions represent Square 
matrices of size 15x15, and the braced expressions represent 
column vectors of size 15x1l. For the reflector, the last 
two terms of equation (64) are zero. Before any operation 
can be performed on equation (64), the nonlinear terms (last 
two terms) must be "linearized" and the term with the 72 
operator must be integrated by parts. 


> 


The nonlinear feedback term, 2n[1 + a Np vy fen] 
was linearized by using predicted values oe 
of the unknown function at time t. These predicted values 
come from the time integration scheme by Franke [7]. 
Basically, the integration scheme utilizes a predictor-cor- 
rector method which predicts values of the unknown function 
at the next time by using the derivatives of the function. 
Adopting the predicted values vy enabled integration over 


Space. The term in equation (64) involving the nonlinear 


feedback can be written as 


Cf NN; N.2n(1 + 1+ Nyy" NoVo' + Pete Ny e¥y5°) Vol] {(c4eu) 5} 
Vol 4 


The last tern of equation (64) describes the delayed 


neutron contribution. In general, 


t 
n y % Ss n-1 \ 
th f 2 )dt' =e Mt, th 1) [e Atha] f et Vs (t"\dt"] 
re) 0 
= n ora 
fe re f ae wan (t!)dt" (65) 

t 
n-1 


Where n is number of time steps, th is the current time, 
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and the] is the previous time. To approximate the integrals 


in equation (65), the predicted values vse were employed 
in a simple trapezoidal rule. The trapezoidal rule was be- 
lieved to be sufficient since very small time steps were uti- 


lized in this work. For example, at time t, , 


(SUM.) = (previous SUM. )e 7° + ee )+ nF) 
j i 2 j ] j 
where 
H = current time step = t. - t, 
-Xt Coy Tet e 
(Sumy) =I" 2 aes va (t!)dt! 
fe) 
= et 
(previous SUM, ) Sg f ane a te et 
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The previous sum is formed by accumulating the trapezoidal 


integration of each time step. In general then, for time th 


- P 
(SUM, ) = (previous SUM. )e + Be - ) 


‘ t 


ne thy (66) 
The delayed neutron term in equation (64) can be expressed 


as 
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In order to bring equation (64) to final form, the V2 
Operator was integrated by parts. Since the flux at the sur- 
face of the reactor is zero, integration by parts yields 


ON. ON. ON. ON. ON. ON. 
‘pee i 


2 es ie ema eee Lis es? 
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Vol Vol 


In vector notation, equation (67) is expressed as 


ON. 
A.J! 
aX 
f aN. oN. oN. ar 
f NV N.dxdydz = - yt oe a See a dxdydz (68) 
Vol Vol 
N. 
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dz 
Using the chain rule produces 
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Letting the 3x3 matrix above be [B'], the V? term becomes 
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Ge - at) 
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where [B']! is the transpose of [B'] . By applying the 
chain rule in equation (47) and using equation (59), it can 


be shown that for this work 


0 
aN aN > ON aN Ns 
1 15 1 15 
Sry %1 s ae X15) or, anaes S X15) 
aa 1 1 
[6°] aN, aN. > ON Ng a Ni 
Ts ean eo aL, Yy5) 
2 
0 3 0 5 h 
[B'y! can be derived from equation (71). 


Note from equation (64) that there are three basic ele- 
ment matrices which are defined as follows after applying 


equation (60): 


Se hitcaeh LS, 
00 


Pe ] j j j j j rad 
(eG, ] 73 Solal OE, 7 aL4) a; = at) ax Le LBs 
-l00 
aN. - 3M. 
(sc - ar) 
1 3 
IN. an. F 
[B'] se 2 a) det[J ]dL,dL,dz 
aN 
a (73) 
te K P P * 
[GGG .] = oP Pes NN, on(14 ae +...4N e046) )det[J dL dL dz 
-100 


(74) 
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Element matrices (4, ;] and [6G 55] are independent of 


time. However, (GGG, 5] is time dependent due to the utili- 
zation of the predicted values vs which changes with time. 


In terms of these three basic element matrices, equation (64) 


is 
[65 ,]{b5) + [GG,,]{cl-v"),} + [G,,]{(c2-v") 5} 
+ [6GG,,] {(c4-v");} + [G,,] {(c5+SuM),} = 0 (75) 


The last two terms of equation (75) are zero for the reflector. 


G. CONSTRUCTION OF THE SYSTEM MATRICES 

The 15x15 coefficient element matrices were calculated 
according to equations (72), (73), and (74); and the results 
were collected element by element into the corresponding sys- 
tem coefficient matrices. The system coefficient matrix 
[BIGG] is developed from [G55] » LBIGGG] from (6G. ] ‘ 
and [BIGH] from [GGG 5] . [BIGG] and [BIGGG] are inde- 
pendent of time and can be constructed once and for all from 
geometry considerations. [BIGH] is dependent on both geom- 
etry and time due to the time dependence of the predicted 
flux utilized in the feedback term. Thus, [BIGH] is recal- 
culated at each time increment. 

Non-zero contributions to a global nodal point I come 
only from adjacent elements sharing that same nodal point I. 
Thus, the system matrices are sparse and banded. The pro- 


cess of assembling contributions from element matrices 
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requires the identification of a local nodal point (i=1,2, 
...,15) with a global nodal point (1=1,2,...NUMNP, where 
NUMNP is the total number of system nodes). This correspond- 
ence between element and global nodes is accomplished via 
the connectivity matrix. 

The formal treatment of the field equations in terms of 


the system coefficient matrices is described by the equation 


[B1GG] {p,} + [BIGGG] {(cl-y),} + [BIGG] {(c2-v),} 


+ [BIGH] {(c4-p) 7} + [BIGG] {(c5+SUM) ,} fn O (76) 


where the system matrices are NUMNP x NUMNP and the column 
vectors are of length NUMNP x 1. However, the direct appli- 
cation of equation (76) requires a large amount of computer 
storage. To take advantage of the sparsity of the system 
matrices, an optimum compacting scheme described by Ref. 8 
was employed. 

The concept behind OCS is simply to store only the non- 
zero terms of a coefficient matrix. OCS requires two integer 
arrays, say JB and NAME, and a vector of non-zero coeffi- 
cients of the square system matrix. For purposes of illus- 
tration, the square system matrix is called B, and the vector 


, ta 


of non-zero coefficients of B is called BB. The i integer 


entry in the NUMNP x 1 JB vector is the number q;- This 


number is defined by 
i=] 
: | by P; s-9=152).25% ,NUMRP (77) 
351 


th 


where P is the number of terms in the i equation. In 
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other words, Pj is the number of nodes that the yuh node 


"sees". JB is therefore a pointer vector of length NUMNP+1 


whose ay term locates the initial position in the BB vector 


of the contributing coefficients to the Kes equation. The 
NUMNP 
NAME vector of length Mx1l, where M = x p 
i=] 
NUMNP successive vector blocks of variable length P. 


a consists of 


> 


integer numbers in the ra block 


j 
of NAME list the P. contributors to the i" equation. The 


Patel... 2tUMNP. The p 


Mx] BB vector contains the real non-zero coefficients of the 
NUMNPXNUMNP B matrix, arranged in the same contiguous block 


th term in the jth 


arrangement as the NAME vector. The j 
block or BB(JB(I) + J-1) is B(I,K), where K = NAME(JB(I)+J-1). 
To illustrate, consider the grid shown in figure 12. The two 


array vectors and the coefficient vector matrix of non-zero 


terms are: 
aa st <1/5, 10, 13, 18, 25, 30, 297 38, aos 
Wee i <1..2.4.512.3,),6.5| --=. 1908.6.5> 
BB = <By 1 9B) 5 5B, 4 5By 5|Bo9 2B 53 Bo] 2B og sBo5|---|Bgq »Bgg »Bg¢ Bg 5> 


In this illustration, NUMNP = 9 and M = 41 [8]. 

In this work, a judicious method of numbering the system 
nodal points was adopted to further reduce computer storage 
requirements. Since the surface boundary nodes of the reac- 
tor represent zero neutron fluxes, the contributions of these 
nodes to interior or non-zero nodes can be discarded. Thus, 
only the interior nodal points need to be considered. These 


non-zero nodes were numbered first in the finite element mesh 
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Figure 12. Sample grid used for illustrating OCS 


57 


Asdmunt stnemels 


296 prisevseur(t 16> bee otto siamese? Sf avughd 


, used so that in the OCS, the number of non-zero nodes (NNZ) 
replaces NUMNP. 

The vectors of non-zero coefficients will be designated 
BIGG, BIGGG and BIGH since the square coefficient matrices 
described in equation (76) were not utilized. To illustrate 
the application of OCS in the system, the following sample 


program is given: 


po 450 I[=1, NNZ 
JBB = JB(1I) 
JE = JB(I+1)-1 
DY(I) = 0.0 
pO 500 J=JBB, JE 
LL = NAME(J) 
OY{1) = BY(I) + BIGG(J)*p(LL) + BIGGG(J)*c1(LL)* 
W(LL) + BIGG(J)*c2(LL)*p(LL) + BIGH(J)*c4(LL)* 
w(LL) + BIGG(J)*c5(LL)*SUM(LL) (78) 
500 CONTINUE 
450 CONTINUE 


where 
Li = nodal point “seen“ by node I 
JE-JBB = total number of nodal points that node I "sees" 
DY(I) = summation of all contributions to node I and 
should sum to zero as stated by the field equa- 
tion 
The assumption of a homogeneous reflector was relaxed in the 
application of equation (78). Interface nodes were assigned 


properties of the core. Therefore, if LL is a reflector 
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node not on the core reflector interface, the last two terms 


of equation (78) are nonexistent. 
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IV. NUMERICAL INTEGRATION 


A. LINE AND AREA INTEGRATION 

The solution of the element matrix equations was achieved 
through numerical integration since an exact closed form solu- 
tion cannot be established. The volume integration was accom- 
plished by using a line integration in the c-direction and an 
area integration in the plane of the triangle. The line inte- 


gral is described by the Gaussian quadrature formula [12] 
] 


n 
ff stcrae sage Daag Pak 6 ty (79) 
“1 as eal 


where nis the number of Gauss integration points 


Hy weighting coefficients 


F(a,) 


the function f(c) evaluated at Gauss point ap 


Table III lists ays Hy and [1] ]. 


The area integration was achieved by the equation 


tak, z 

J m m m 
ia f(L,,L,,L,)dL,dLy = 2) wif(L, sty »L3 ) (80) 
> 58 , 


where m is the number of area integration points and Wo 
are the weights. The numerical integration points for the 


area integration are given in Table IV which was extracted 
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TABLE {11 


Abscissae and Weight Coefficients 
of the Gaussian Quadrature Formula 


] 
ff tleyde = ED Atay) 
= k=1 
+a H 
he 72 
0:°57735 02691 89626 1-00000 00000 00000 
: n= 3 
0:77459 66692 41483 0°55555 55555 55556 
0:00000 00000 00000 0-88888 88888 88889 
n=4 
0-86113 63115 94053 0-34785 48451 37454 
0-33998 10435 84856 0:65214 51548 62546 
n=5 
0:90617 98459 38664 0:23692 68850 56189 
0:53846 93101 05683 0-47862 86704 99366 
0:00000 00000 00000 0:56888 88888 88889 
n=6 
0:93246 95142 03152 0:17132 44923 79170 
0:66120 93864 66265 0-36076 15730 48139 
0:23861 91860 83197 0°46791 39345 72691 
n= 7 
0:94910 79123 42759 0:12948 49661 68870 
0:74153 11855 99394 0:27970 53914 89277 
0-40584 51513 77397 0-38183 00505 05119 
0-00000 00000 00000 0:41795 91836 73469. 
n=8 
0-96028 98564 97536 0-:10122 85362 90376 
0:79666 64774 13627 0:22238 10344 53374 
0°52553 24099 16329 9:31370 66458 77887 
0:18343 46424 95650 0:36268 37833 78362 
n=9 
0:96816 02395 07626 0:08127 43883 61574 
0:83603 11073 26636 0-18064 81606 ©4857 
0:61337 14327 00590 0:26061 06964 v2935 
0-32425 34234 03809 0:31234 70770 40002 
0:00000 00000 00000 0:33023 93550 01260 
n= 10 
0:97390 65285 17172 0:06667 13443 08688 
086506 33666 88985 014945 13491 50581 
0:67940 95682 99024 0:21908 63625 15982 
0:43339 53941 29247 0:26926 67193 09996 
0:14887 43389 81631 0:29552 42247 14753 
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TABLE IV 


Numerical Formulas for Triangles 


Order Fig. Error Points Triangular Weights 
Co-ordinates 2W, 
Linear ADS R = 0(h?) a A kert 1 
a 4,4,0 4 
Quadratic R = 0(h?) b 0, $.4 4 
ce =. 4,0, } 5 
EN a 4.44 -¥ 
Cubic Je Q R = Oh4) b i fa FH 
IN dts aH 
a AS 25 
b 4, 4.0 
QW, G 0.4.4 \ 20 
Cubic d R = \{h*) d 40.4 
Se Se py, 
} ek 
g 0,0.1 
ao nat 0-225 
) b 4. BB 
ce: «Bland, 0:13239415 
Quintic R = QA") d B,. By. 2, 
A> e 22, Bz, B: 
<—— | eee ee 0-12593918 
g 
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B;. Bz, 22 


with 

2, = 0:08971587 
B, = 047014206 
x, = 0:79742699 
B, = 0:10128651 


from Ref. 11. The volume integration of the function 


F(L,,L,5L352) using numerical integration is therefore 


17 I-L; 
Sf f ee atgetg eran yat ac ~ 
aq? toro | 
n m 
m m m ok 
Re Wy 2 minh Ly tho sky seu (81) 
where tf = Hy . In the manner of equation (81), the element 
matrices become 
h n m m m m =k os 
[65,] = 5 2» Wy ay We N5(L, »Ly sb3 5% )N,det[J Pe (82) 
h n m m m m 
(ac. ] ez 2 Wy x Wo F (Ly »Ly sla 5% \ netted. 21, 033) 
= m= 
h n m m m ma - 
[occ] = > 5 e Wo WoT Ah by bse pedecle 2oaier 
k=] m=] 


wmyere F “and T are easily derived from equations (73) and 
(74), respectively. Note that N. and det[J ] are also 
evaluated at each integration point. They are not shown as 
Such merely for the sake of convenience in writing the equa- 


EONS . 


B. NUMBER OF INTEGRATION POINTS 

It is difficult to estimate the number of integration 
points required for good accuracy due to the complexity of 
the functions involved. The basic rule that the best number 
of integration points is found by trial and experience was 


adopted. For the [G5,] element matrix, the numbers involved 
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can be approximated. This was done by letting the determi- 
nant of the Jacobian equal to twice the area of the curved 
triangle (checking det[J*] at each integration point showed 
that this assumption was not too unreasonable). With det[J*] 
outside the integration process, [G55] can be solved in 
closed form by integrating out ¢ from -1] to +1 and then apply- 
ing the closed form equation [12] 
m m m,!m,!m,! 
2 L, 3-d(Area) = 2A Tate sa FOIE 

(85) 
where mj> M5, Mz are positive integer exponents and A is 
the area of the triangle. 

Five test points within the 15x15 element matrix were 
selected, as listed in Table V. The values obtained from 
the application of equation (85) are also given in Table V. 
Three sets of area integration points, each with a different 
number of c Gauss points, were used. These three sets of 
area integration points, given in Table IV, are: 

me cuore order (4 points) 

2. c¢ubte order (7 points) 

we guintic order (7 points) 

Using each of the three integration points above with differ- 
ent c Gauss points in equation (82) produced the results 
obtained in Table V. From these results, the quintic order 
area integration points were selected for the element matrix 
[G55] with three c Gauss points in the prism axes. This 


set of integration points was also used for [GGG 55]. 
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TABLE V. 


Selection of Integration Points for [G5] 


Area 
cS points oints G(1,] G(2,2 G(1,9 G(6,9 G(9,9 


ee cubic 1415 5278 -2756 5072 8536 
(4 pts) 

5 z 1817 5278 -3170 5072 10240 

"i ‘ Isl? 5278 -3170 5072 10240 

3 cubic 3248 8247 -3114 4960 10243 
(7 pts) 

2 , 3249 8247 -3114 4960 10240 

7 3249 8247 -3114 4960 10240 

2 quintic 2464 6600 -3144 5020 10240 
(7 pts) 

5 si 2464 6600 -3144 5020 10240 

7 : 2464 6600 -3144 5020 10240 

Approximated values 2500 6700 -3142 5026 10053 
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The numbers for the (GG, 5] element matrix could not 
be approximated. The best that could be done was to obtain 
an idea of the order of magnitude of this element matrix. 
Towards this end, a linear triangular element was assumed. 
Using linear approximation, the order of magnitude was found 
to be about hah Using the three different area integration 
points mentioned above with varying 7% Gauss points in equa- 
tion (83) yielded the results given in Table VI. Further 
checks on the [665 ,] element matrix showed that the cubic 
order with four area integration points yielded the desired 
order of magnitude. Thus, the fourth order cubic with five 
¢ Gauss points was employed for the (GG, 5] element matrix. 
A note should be mentioned here in regards to the vast dif- 
ference in results obtained for (6G, 5] using different 
area integration points. Most likely, it was due to the 
FB"] and fan! matrices which required the inversion of 
ar, ; aL, , ete. Or perhaps 1t was caused by the nature of 


the hybrid element used in this work. In any case, further 


investigation is warranted in this area. 
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TABLE VI 


Selection of Integration Points for [GG 5] 


Area 
oints oints 


cubic 
(4 pts) 


quintic 
(7 pts) 


GG(1,9) 
18.4 


22.0 
eeu 


-1.03x10* 


~1,03x10" 
-1.03x10° 


-52.6 


-52.6 
-52.6 


GG(3,3 
170.1 


1776 
1776 
2. PIxro 


Zt into 


2.11x10" 
1.27x104 


1.27x107 
1.27x10° 
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GG(14,9 
-165 


-186 


-186 
1.05x10/ 


105x107 


1,05x10/ 
6 64x10" 
~6.64x10° 


-6.64x10" 


GG(8,8 
165.9 


195.9 


195.9 
8. 42x10" 
8.42x10/ 


8.41x10" 


- 


6.84x10 


6.84x10° 
6.84x10" 


GG(15,15 
106.1 


106.1 
106.1 
2056 


2056 


2056 
1.51x10° 
1.51x10° 


“1.51x10° 


V. TEST PROBLEMS AND RESULTS 


The reactor was subjected to uniform and local perturba- 


tions in the form of a ramp input described by 
m)= 5, + at (86) 


where a is the change in Le per unit time and oo is 
the critical fission cross-section. ae must first be ob- 
tained before the perturbations can be applied. This was ac- 
complished by trial and error until a stationary solution 
was reached. For mesh l, xo was found to be 0.0057360 
per cm. No attempt was made to find zt for mesh II due 
to time limitations. As such, the test problems outlined be- 
low were applied to mesh I. Future work is planned to apply 
the test problems on mesh II. 

The following perturbations were applied: 

a) Uniform perturbation of 10 dollar of reactivity per 


second: 
(r,t) = ia + ees, in the core 
where a = 0.005893/cm-sec 


b) Local perturbation at the core center of 100 dollar 


of reactivity per second: 


* a 
ryt)c= zg +até(r) ; *dn® thec tore 
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where r, is (0,0,0) and a = 0.015407/cm-sec 
c) Local, off-center perturbation: 


* = 
Ze(r,t) a a,to(r,) > 43 = 1,228. i” the care 


where 
ry = (0,60,40) 
a, = 0.015407/cm-sec = 100 dollar per second 
do = 0.008123/cm-sec = 50 dollar per second 
a3 = 0.005894/cm-sec = 10 dollar per second 


emree est points, (0,0,0), (60,0,0),./and (-60.,0,80). 
were selected to trace the neutron time history. For cases 
a) and b), the neutron flux was plotted at each test point 
during transience. This 1s shown in figures 13 and 14. 
Case c) involved three ramp inputs and were conducted for 
both the linear and nonlinear reactor equations. The linear 
and nonlinear responses were compared at each test point for 
each ramp input and are illustrated in figures 15 through 23. 
The radial and axial flux distributions at time t = 0.0123 
second were plotted for the steady state and the 100 dollar 
perturbation case. These are embodied in figures 24 and 25. 
Finally, a neutron flux early time history between mesh I 
and mesh II was plotted to check the effect of using a finer 
element mesh. The result is portrayed in figure 26. 

Figures 13 and 14 revealed a clear space dependence of 
the neutron flux during transience as expected. Figures 15 


through 23 demonstrated the effects of temperature feedback 
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Figure 13. Neutron flux transient history at three 
test points with a uniform perturbation 
of 10 dollar of reactivity per second. 
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Figure 14. Neutron flux transient history at various 
test points for a local central perturbation 
ef “100 “dol Tar sec “of “regcetivrey. 
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neutron flux 


Figure 15. 


time (sec) 


A - nonlinear neutron response 
B - linear neutron response 


Linear and nonlinear fluxes at (0,0,0) due 
to a local perturbation of 100 dollar/sec of 
reactivity at (0,60,40). 


72 


(Sez) smd 


sen 1 noytwsn Vsentinen - A 
Non2st puysusp vesu?l - a 


a ctariek J6 250? ‘sent toon bat yeant 
Yo rfeb OOF Fo Ag nal tay [adolf § oF 
-Co-0) Js ysrvreoeey 


2 


10 
B 
11 
= 16 A 
i 
= 
° 
Ls 
+ 
=| 
<8) 
iS 
10! 
10° 
ror 10=" ro" 10" 


time (sec) 


A - nonlinear neutron response 
B - linear neutron response 


Figure 16. Linear and nonlinear fluxes at (60,0,0) due 
to a local perturbation of 100 dollar/sec at 
(0560240): 
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A - nonlinear neutron response 
B - linear neutron response 


Linear and nonlinear fluxes at (-60,0,80) 
due to a local perturbation of 100 dollar/sec 
at (0,60,40). 
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Figure 18. 


time (sec) 


A - nonlinear neutron response 
B - linear neutron response 


Linear and nonlinear fluxes at (0,0,0) due 
to a 50 dollar/sec local perturbation at 
(0,60,40). 
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Figure 19. 
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time (sec) 
A - nonlinear neutron response 
B - linear neutron response 


Linear and nonlinear fluxes at (60,0,0) 
due to a 50 dollar/sec local perturbation 
at (0,60,40). 
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Figure 20. 
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A - nonlinear neutron response 
B - linear neutron response 


Linear and nonlinear fluxes at (-60,0,80) 
due to a 50 dollar/sec local perturbation 
at, (0.60.40). 
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A - nonlinear neutron response 
B - linear neutron response 


Linear and nonlinear fluxes at (0,0,0) due 
to a 10 dollar/sec local perturbation at 
(0,60,40). 
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Figure 22. Linear and nonlinear fluxes at (60,0,0) due 
to a 10 dollar/sec local perturbation at 
(0,60,40). 
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Figure 23. Linear and nonlinear fluxes at (-60,0,80) 
due to a 10 dollar/sec local perturbation 
at {0.6040;). 
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Figure 24. Radial flux distribution for the steady state 
and 100 dollar/sec local perturbation. 
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Figure 25. Axial flux distribution for the steady state and 
100 dollar/sec local perturbation. 
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ragure 26. Early time history of the neutron Flux.at core 
center using mesh I and mesh II. 
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and delayed neutrons on the flux. To obtain an idea of the 
difference in magnitude between the linear and nonlinear 

feux for the 100 dollar local perturbation, at time t = 107? 
second, the linear case predicted a neutron flux at the core 


center of 2.024 x 10!@ neutrons per cm? per sec. The non- 


linear reactor predicted 1.85 x 19/9 neutrons per cm? per 
sec. Although the numbers are not sufficiently accurate due 
to the crudeness of the finite element mesh, the significant 
difference in the order of magnitude between the linear and 
nonlinear neutron flux leads to the belief that the three- 
dimensional finite element model utilized is predicting the 
correct trend of neutron flux behavior. 

The radial flux distribution for the steady state shown 
in figure 24 portrays the expected flux distribution. The 
radial flux distributions for the local perturbation case 
show a skew distribution at the point of perturbation. The 
axial flux distribution at steady state is very much symme- 
tric about the center. Again, the expected skew at the 
point of perturbation is there, as shown in figure 25. 
Figure 26 gives an indication that the ah for finer 
meshes is higher. However, curve B of figure 26 should be 
extended to longer times to verify this hypothesis. Curve B, 
as plotted in figure 26, used four hours of computer time 


employing the H-compiler of the IBM 360/67. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


The finite element mesh employed here is crude, and, 
thus, results that were obtained should be considered as posi- 
tive indicators rather than numerically conclusive facts. The 
trend of neutron flux behavior is the major thrust of this 
work. From the results obtained, it is concluded that the 
expected patterns of neutron behavior as predicted by the 
three-dimensional finite element model used do occur. These 
patterns were best demonstrated by the differences between 
the linear and nonlinear flux responses and by the spatial 
flux distribution at steady state and during local perturba- 
tion. The three-dimensional quadratic finite element model 
utilized should produce better results by resorting to a 
finer mesh. The draw-back to a finer mesh, of course, is 
the significant increase in computer time and storage re- 
quirements. Once accurate results are obtained through 
finer element meshes, comparisons between three- and two- 
dimensional models can be attempted. 

A mesh generator was not developed for this problem. It 
is recommended that this be done to ease the transition from 
One mesh to another and to minimize human error. In addi- 
tion, a similar calculation using a three-dimensional linear 
element should be performed to corroborate the results ob- 
tained here. It should be particularly noted that this type 


of problem is highly sensitive to the fission cross-section. 
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In the search for Xe » it was necessary to adjust 2 to 


Fo 
the sixth decimal place. There is no exact method of deriv- 
ing i due to the highly nonlinear aspect of the problem; 
therefore, trial and error must be used. Finally, the 
Gaussian quadrature used to determine the element matrices 
should be investigated. The cause of the differences in 
values obtained by using different number of integration . 
points should be established. Was it due to the integration 
points, the coordinate transformation, or the element shape 
functions? This is an important question upon which future 


numerical results could be based. 
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APPENDIX A 


MESH I CONNECTIVITY AND COORDINATES 


The connectivity matrix for the 128-element mesh is as follows 


1l 12 13 14 15 
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The coordinates of the 128-element mesh are as follows 
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APPENDIX B 


MESH II CONNECTIVITY AND COORDINATES 
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The connectivity matrix for the 192-element mesh 
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The coordinates of the 192-element mesh are as follows 
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